Numerical experiments of adjusted BSSN systems for controlling constraint violations 
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We present our numerical comparisons between the BSSN formulation widely used in numerical 
relativity today and its adjusted versions using constraints. We performed three testbeds: gauge- 
wave, linear wave, and Gowdy-wave tests, proposed by the Mexico workshop on the formulation 
problem of the Einstein equations. We tried three kinds of adjustments, which were previously 
proposed from the analysis of the constraint propagation equations, and investigated how they 
improve the accuracy and stability of evolutions. We observed that the signature of the proposed 
Lagrange multipliers are always right and the adjustments improve the convergence and stability of 
the simulations. When the original BSSN system already shows satisfactory good evolutions (e.g., 
linear wave test), the adjusted versions also coincide with those evolutions; while in some cases (e.g., 
gauge-wave or Gowdy-wave tests) the simulations using the adjusted systems last 10 times as long 
as those using the original BSSN equations. Our demonstrations imply a potential to construct a 
robust evolution system against constraint violations even in highly dynamical situations. 

PACS numbers: 04.20.-q, 04.25. Dm, 04.25.-g 



I. INTRODUCTION 

Numerical integration of the Einstein equations is the 
only way to investigate highly dynamical and nonlinear 
gravitational space-time. The detection of gravitational 
wave requires templates of waveform, among them merg- 
ers of compact objects are the most plausible astrophysi- 
cal sources. Numerical relativity has been developed with 
this purpose over decades. 

For neutron star (NS) binaries, a number of scien- 
tific numerical simulations have been done so far, and 
we are now at the level of discussing the actual physics 
of the phenomena, including the effects of the equa- 
tions of state, hydrodynamics, andgeneral relativity by 
evolving various initial data [H, UHi H HI- Mergers 
of black holes (BHs) are also available after the break- 
through by Pretorius [6j in 2004. Pretorius's implemen- 
tation had many novel features in his code; among them 
he discretizes the four-dimensional Einstein equations di- 
rectly, which is not a conventional approach so far. How- 
ever, after the announcements of successful binary BH 
mergers by Campanelli et al. [3] and Baker et al. @ 
based on the standard 3+1 decomposition of the Einstein 
equations, many groups be gan producing interesting re- 
sults [I EE El El EE EllS ES E3, UM- The merger 
of NS-BH binary simulations has also been reported re- 
cently, e.g. 

Almost all the groups which apply the above conven- 
tional approach use the so-called BSSN variables together 
with "1 + log" -type slicing conditions for the lapse func- 
tion and 'T-driver" type slicing conditions for the shift 
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function. BSSN stands for Baumgarte-Shapiro [2fJ and 
Shibata-Nakamura [21| . the modified Arnowitt-Deser- 
Misner formulation initially proposed by Nakamura [22] • 
(The details are described in i jll Al ) There have al- 
ready been several efforts to explain why the combi- 
nation of this recipe works from the point of view of 
the well- pos edness of the partial differential equations 
(e.g. [23L l24j). However, the question remains whether 
there exists an alternative evolution system that enables 
more long-term stable and accurate simulations. The 
search for a better set of equations for numerical inte- 
grations is called the formulation problem for numerical 
relativity, of which earlier stages are reviewed by one of 
the authors [251 ]. 

In this article, we report our numerical tests of mod- 
ified versions of the BSSN system, the adjusted BSSN 
systems, proposed by Yoneda and Shinkai [26] . The idea 
of their modifications is to add constraints to the evolu- 
tion equations like Lagrange multipliers and to construct 
a robust evolution system which evolves to the constraint 
surface as the attractor. Their proposals are based on the 
eigenvalue analysis of the constraint propagation equa- 
tions (the evolution equations of the constraints) on the 
perturbed metric. For the ADM formulation, they ex- 
plain why the standard ADM does not work for long-term 
simulations by showing the existence of the constraint vi- 
olating mode in perturbed Schwarzschild space-time [27| . 
For the BSSN formulation, they analyzed the eigenval- 
ues of the constraint propagation equations only on flat 
space-time [261 ] . but one of their proposed adjustments 
was immediately tested by Yo et al. [28j for the numeri- 
cal evolution of Kerr-Schild space-time and confirmed to 
work as expected. (The details arc described in SQlBj) 

Our numerical examples are taken from the proposed 
problems for testing the formulations of the Mexico 
Numerical Relativity Workshop 2001 participants [29| . 
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which are sometimes called the Apples- with- Apples test. 
To concentrate the comparisons on the formulation prob- 
lem, the templated problems are settled so as not to 
require technical complications; e.g., periodic boundary 
conditions are used and the slicing conditions do not re- 
quire solving elliptical equations. Several groups already 
reported their code tests using these Apples tests (e.g. 
[3Q , I3II , H2| ) , and we are also able to compare our results 
with theirs. 

This article is organized as follows. We describe the 
BSSN equations and the adjusted BSSN equations in 
Sec. Ill Al and III Bl We give our three numerical test 
problems in Sec. IIIII Comments on our coding stuff are 
in Sec. [IVj Sec. |V] is devoted to showing numerical re- 
sults for each testbeds, and we summarize the results in 
Sec. ED 



widely used among numerical relativists. 

The idea of the BSSN formulation is to introduce aux- 
iliary variables to those of the Arnowitt-Deser-Misner 
(ADM) formulation for obtaining longer stable numerical 
simulations. The basic variables of the BSSN formulation 
are (0, 7^ , K, Aij , T*), which are defined by 
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II. BASIC EQUATIONS 

A. BSSN equations 

We start by presenting the standard BSSN formula- 
tion, where we follow the notations of [20(, which are 



where (7^ , Kij ) are the intrinsic and extrinsic ADM 3- 
metric. The conformal factor <p is introduced so as to 
set 7 = det[7y] as unity, Aij is supposed to be trace- 
less, and r l is treated independently in evolution equa- 
tions. Therefore these three requirements turn into the 
new constraints [below (|2T6|l - (|238|) ]. 

The set of the BSSN evolution equations are 
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where Di is the covariant derivative with respect to the 3- metric 7y and TF means trace- free operation, i.e., Hj F 
Hij — ^jijH k k . The Ricci tensor is computed with the conformal connection T l as 



R% = -2A£»j0 - 2^D k D k 4> + 4D l( j>D^ - A%D k 4>D k ^ 
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where Di is a covariant derivative associated with 7^ . the Hamiltonian and momentum constraint equations, 
Similarly to the ADM formulation, this system has are expressed in terms of the BSSN basic variables and 
constraint equations. The two "kinematic" constraints, are written as 
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Additionally, the BSSN formulation requires three "alge- 
braic" constraint relations; 
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where f|2 . and (|2 . 1 7[) are from the definitions of 
and (|2.4|) . respectively. Equation (|2.18|) is from the re- 
quirement on 7. 

These five constraints are, theoretically, supposed to be 
zero at all times; therefore they can be used to modify the 
dynamical equations. For example, Alcubierre et al. [34| 
announced that the replacement of the terms in (|2.10j) 
using the momentum constraint drastically changes the 
stability feature. Actually, such replacements of terms 
using constraints are applied (with/without intentions) 
in many terms in (|2.6D - (|2.10[) . which are expressed as 
Eqs. (2.27)-(2.31) in j2j|. 

Alcubierre et al. [35| also pointed out that the re- 
definition of Aij by 



A, 



A,. 



trA 



(2.19) 



during the time evolutions improves the numerical sta- 
bility. This technique again can be understood as the 
trace-out of the .A-constraint (|2.17|) from the evolution 
equations. In our numerical code, we do not apply this 
technique because we recognize the trace-free property 
as the new constraint A in the BSSN system, and our 
purpose is to construct a system preventing the violation 
of constraints. 

Recently, several groups applied artificial dissipation 
(e.g. [36] ]) to obtain stable evolutions (see, e.g. [HI, [H, 
137]). We, however, do not introduce such dissipations 
in our code, since we try to clarify the difference of sta- 
bility from the viewpoint of formulations of the Einstein 
equations. 



B. Adjusted BSSN systems 

To understand the stability property of the BSSN sys- 
tem, Yoneda and Shinkai [2(| studied the structure of the 
evolution equations, (|2.6p - (|2.10p , in detail, especially how 
the modifications using the constraints, (|2.14p - (|2.18p . af- 
fect to the stability. They investigated the signature of 
the eigenvalues of the constraint propagation equations 
(dynamical equations of constraints), and explained that 



the standard BSSN dynamical equations are balanced 
from the viewpoints of constrained propagations, includ- 
ing a clarification of the effect of the replacement using 
the momentum constraint equation. 

Moreover, they predicted that several combinations 
of modifications have a constraint-damping nature, and 
named them adjusted BSSN systems. (Their predictions 
are based on the signature of eigenvalues of the constraint 
propagations, and the negative signature implies a dy- 
namical system which evolves toward the constraint sur- 
face as the attractor.) 

Among them, in this work, we test the following three 
adjustments: 

1. An adjustment of the A-equation with the momen- 
tum constraint: 



d t A 



K A aD(iMj), 



(2.20) 



where ka is predicted (from the eigenvalue analy- 
sis) to be positive in order to damp the constraint 
violations. 

2. An adjustment of the 7-equation with Q constraint: 
daij = + ^a% (i D j) g k , (2.21) 



with k~ < 0. 



3. An adjustment of the T-equation with Q constraint: 
d t f l = dft i + K t aQ i . (2.22) 



with Kp < 0. 



These three adjustments are chosen as samples of "best 
candidates", Eq. (4.9)-(4.11) in 26]. The term "best" 
comes from their conjecture on the eigenvalue analysis 
of the constraint propagation matrix; that is, (a) all the 
resultant eigenvalues from above adjustments can be less 
than or at most equal to zero, which indicates the de- 
cay of constraint errors, and (b) the resultant constraint 
propagation matrix is diagonalizable, which guarantees 
the predictions of above eigenvalue analysis (see Table II 
in [26j). However, since above eigenvalues include zero 
elements and also above analysis assumes a linearly per- 
turbed metric about the flat space-time, the effects of 
the adjustments (|2.20p - (|2.22p need to be demonstrated 
via numerical experiments. 
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III. NUMERICAL TESTBED MODELS 



where 



Following the proposals of the Mexico Numerical Rela- 
tivity Workshop [29j, we perform three kinds of tests. In 
this section, we explicitly give some details of the models. 



A. Gauge-wave testbed 

The first test is the trivial Minkowski space-time, but 
sliced with the time-dependent 3-metric, which is called 
the gauge- wave test. The 4- metric is obtained by coordi- 
nate transformation from the Minkowski metric as 



ds 2 



-Hdt 2 + Hdx 2 + dy 2 + dz 2 , 



where 



H = H(x -t) = l- Asin 



2tt(x - t) 



(3.1) 



(3.2) 



which describes a sinusoidal gauge wave of amplitude A 
propagating along the x-axis. The non-trivial extrinsic 
curvature is 



K, 



irA 



1 + A sin 



2ir(x-t) 



(3.3) 



Following [29l | . we chose numerical domain and parame- 
ters as follows: 

• Gauge- wave parameters: d — 1 and A = 10~ 2 

• Simulation domain: 0.5, 0.5], y = z = 

• Grid: x l = —0.5 + (n — ^)dx with n = 1, • • • 50p, 
where dx = l/(50p) with p = 2,4, 8 

• Time step: dt = 0.25dx 

• Boundary conditions: Periodic boundary condition 
in x direction and planar symmetry in y and z di- 
rections 



• Gauge conditions: 

d t a = -a 2 K, ft = 0. 



(3.4) 



The ID simulation is carried out for a T — 1000 crossing- 
time or until the code crashes, where one crossing-time 
is defined by the length of the simulation domain. 



B. Linear wave testbed 

The second test is to check the ability of handling a 
travelling gravitational wave. The initial 3-metric and 
extrinsic curvature Kij are given by a diagonal pertur- 
bation with component 



b = A sin 



2tt(x - t) 



(3.6) 



for a linearized plane wave traveling in the ir-direction. 
Here d is the linear size of the propagation domain and A 
is the amplitude of the wave. The non-trivial components 
of extrinsic curvature arc then 



1 



1 



K yy = ~7;dtb, K zz = -d t b. 



(3.7) 



2 2 
Following [2^|, we chose the following parameters: 

• Linear wave parameters: d = 1 and A = 10~ 8 

• Simulation domain: xg[~ 0.5, 0.5], y — 0, z = 

• Grid: x l = —0.5 + (n — ^)dx with n = 1, • • • 50p, 
where dx = l/(50p) with p = 2, 4, 8 

• Time step: dt = 0.25dx 

• Boundary conditions: Periodic boundary condition 
in x direction and planar symmetry in y and z di- 
rections 

• Gauge conditions: a = 1 and (3 l — 

The ID simulation is carried out for a T = 1000 crossing- 
time or until the code crashes. 



C. Collapsing polarized Gowdy-wave testbed 

The third test is to check the formulation in a strong 
field context using the polarized Gowdy metric, which is 
written as 

ds 2 = t- l ' 2 e x ' 2 {-dt 2 + dz 2 ) + t(e p dx 2 + e- p dy 2 )(3.8) 

Here time coordinate t is chosen such that time increases 
as the universe expands. Simple forms of the solutions, 
P and A, are given by 

P = J (27rf)cos(27rz), (3.9) 
A = -27rfJ (27rf)Ji(27rf)cos 2 (27rz) 

+ 27T 2 t 2 [J 2 (27rf) + J 2 (27rf)] 

1, 



^ (27r) 2 [J 2 (27r) + J 2 (2^)] 
-2^(2^)^(2^)], 



(3.10) 



where J n is the Bessel function. The non-trivial extrinsic 
curvatures are then 



ds 2 



-dt 2 + dx 2 + (1 + b)dy 2 + (1 - b)dz 2 , (3.5) 



K 



K 7 . 



_l t V4 e -A/4 e P (1+tP)t)) (3 . n) 

_ItV4 e - A /4 e --P(i_tP il ) i (3.12) 

It-V4 e A/4 (t -l_ At) (313) 
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According to [29j , the new time coordinate t, which satis- 
fies harmonic condition, is obtained by coordinate trans- 
formation as 

t(r) = ke CT , (3.14) 

where c and k are arbitrary constants. Using this free- 
dom, we can set the lapse function in the new time co- 
ordinate to be unity at the initial time. Concretely, we 
set 

to = r ~ 9.8753205829098, 

c - 0.0021195119214617, (3.15) 
k ~ 9.6707698127638, 

where to is the initial time. Following [29l |. we perform 
our evolution in the collapsing (i.e. backward in time) 
direction. Parameters are chosen as follows: 

• Simulation domain: z £ [—0.5, 0.5], x = y = 

• Grid: z — —0.5 + (n — ^)dz with n = 1, • • • 50p, 
where dz = l/(50p) with p = 2,4,8 

• Time step: dt = 0.25dz 

• Boundary conditions: Periodic boundary condition 
in z-direction and plane symmetry in x- and y- 
directions 

• Gauge conditions: the harmonic slicing (|3.4j) and 

/& = 

The ID simulation is carried out for a T — 1000 crossing- 
time or until the code crashes. 



IV. THE CODE 
A. Code description 

We have developed a new numerical code based 
on the adjusted BSSN systems. The variables are 
(</>, Jij,K, Aij, T l ), and the evolution equations are (|2.6p ~ 
(f2~TUj) with/without adjustment (f2T20|) . (f2~2Tj) , and/or 
(|2.22[) . The time-integration is under the free- evolution 
scheme, and we monitor five constraints, (|2.14|) - (|2.18|) , to 
check the accuracy and stability of the evolutions. 

Our time-integration scheme is the three-step iterative 
Crank-Nicholson method with centered finite difference 
in space [39l | . This scheme should have second-order con- 
vergence both in space and time, and we checked its con- 
vergence in all the testbeds. 

As we have already mentioned in the end of §11 A, we 
do not apply the trace-out technique of , (|2.19p in our 
code. 

We also remark on our treatment of the conformal con- 
nection variable f 1 . As was pointed out in [H| , it is better 
not to use T l in all the evolution equations. We surmise 
this is because the amplification of the error due to the 



discrepancy of the definition (|2.5p . i.e., the accumula- 
tions of the violations of ^'-constraint (|2.16[) . Therefore, 
we used the evolved T 1 only for the terms in (|2.10|) and 
(|2.13[) . and not for other terms, so as not to implicitly 
apply the ^'-constraint in time evolutions. 



B. Debugging procedures 

It is crucial that our code can produce accurate re- 
sults, because the adjustment methods are based on 
the assumption that the code represents the BSSN sys- 
tem (|2.6p - (|2.10p accurately. We verified our code by com- 
paring our numerical data with analytic solutions from 
the gauge-wave and Gowdy-wave testbeds in Sec. IIIII 
The actual procedures are as follows: 

1. Evolve only one component, e.g. A xx , numerically, 
and express all the other components with those of 
the analytic solution. In this situation, the origin 
of the error is from the finite differencing of the 
analytic solution in the spatial direction and from 
that of the numerically evolved component (A xx ) 
both in spatial and time directions. We checked 
the code by monitoring the difference between the 
numerically evolved component (A xx ) and its ana- 
lytic expression. This procedure was applied to all 
the components one by one. 

2. Evolve only several components, e.g., A xx and T x , 
numerically, and express the other components by 
the analytic solution. The error can be checked by 
a procedure similar to the one above. 

3. Evolve all the components numerically, and check 
the error with the analytic solution. 

We repeated these procedure three times by switch- 
ing the propagation directions (x, y, and z-directions) 
of gauge-wave and Gowdy-wave solutions. We also ap- 
plied these procedures in a 2D test [29[ , and checked the 
off-diagonal component. 



C. Error evaluation methods 

It should be emphasized that the adjustment effect has 
two meanings, improvement of stability and of accuracy. 
Even if a simulation is stable, it does not imply that the 
result is accurate. We judge the stability of the evolution 
by monitoring the L2 norm of each constraint, 



\\SCMt)= /-£(C(t; W )) 2 , (4.1) 

where N is the total number of grid points, while 
we judge the accuracy by the difference of the met- 
ric components gij (t; x, y, z) from the exact solution 



6 



7iJ xact ' ) (i; x, y, z), higher resolution cases. However, it is also true that all 

errors are still under the errors of the plain BSSN system. 
/ i / (exact) \ 2 Therefore, we conclude that this adjusted system shows 

\\^9ij\\2{t) = . — \ 9ij ~ 9ij J ■ (4-2) a weaker instability than the plain system. 



D. Magnitude of k 

Adjusted systems, (|2.20l) - (|2.22l) . require to specify the 
parameter n. From the analytical prediction in [261 ] we 
know the signature of k, but not for its magnitude. By 
definition of the adjustment terms in Eq. (|2.20[) - (|2.22[) . 
applying small k should produce the close results with 
those of the plain system. On the contrary, the large 
k system will violate the Courant-Friedrich-Lewy condi- 
tion (4^|. Hence, there exists a suitable region in the 
adjustment parameters. 

At this moment, we have to chose k experimentally, 
by observing the life-time of simulations. The value of /?, 
used in our demonstrations, is one of the choices of which 
the adjustment works effectively in all the resolutions. 

V. NUMERICAL RESULTS 
A. Gauge-wave test 

1. The plain BSSN system 

As the first test, we show the plain BSSN evolution 
(that is, no adjustments) in Fig. [T] for the gauge-wave 
test. In Fig. [TJ the L2 norms of the Hamiltonian and 
momentum constraints (|4. 1[) are plotted as a function of 
the crossing-time. The second-order convergent nature 
is lost at an early time, the 20 crossing-time, and the 
simulation crashes at about the 100 crossing-time. The 
poor performance of the plain BSSN system for the gauge 
wave test has been reported in [3l[ (see their Fig. 8). This 
drawback, on the other hand, can be overcome if one uses 
the fourth-order finite differencing scheme, an example of 
which can be seen in [33] (see their Fig. 2). 

2. Adjusted BSSN with A-equation 

We found that the simulation lasts 10 times longer 
with the adjustment in the A-equation using the mo- 
mentum constraint (|2.20[) . Figure [3] shows the L2 norms 
of the Hamiltonian and Momentum constraints in the 
same style as in Fig. [T] The adjustment parameter is set 
at ka = 0.005 for this plot. We obtain almost prefect 
overlap of the rescaled Hamiltonian constraint for 200 
crossing-times and almost perfect overlap in the momen- 
tum constraint for 50 crossing-times; there apparently 
improve the results of the plain BSSN system (see Fig.[T]). 
We show the plots until the 1000 crossing-time, there we 
observe the growth of the error both in later time and in 



3. Adjusted BSSN with T -equation 

The case of the adjustment of the T-equation using the 
Q constraint (|2.22[) is shown in Fig. [3] 

The adjustment parameter is set at Kf — —0.1. We 
find that the second-order convergence breaks down near 
the 40 crossing-time under the Momentum constraint, 
which is almost the same as with the plain BSSN system. 
However, the convergence of the Hamiltonian constraint 
is improved, i.e., it continues to the near 55 crossing-time. 
The life-time of the simulation is almost the same as that 
of the plain BSSN system. 



4- Adjusted BSSN with ^/-equation 

We also tested the cases of the adjustment of the 7- 
equation using the Q constraint, (|2.21[) . We again ob- 
served the effects of the adjustment on its stability and 
accuracy but found a rather small effect compared to the 
cases of the adjustments of (|2.20[) or (|2.22p . up to our 
trials of the parameter range of k 7 . Therefore we omit 
showing the results. 



5. Evaluation of Accuracy 

For evaluating the accuracy, we prepare Fig. Ufa), in 
which we plot the L2 norm of the error in j xx , (|4.2p , with 
the function of time. Three lines correspond to the result 
of the plain BSSN system, A-eq. adjusted, and f-eq. 
adjusted BSSN system, respectively. The T-adjustment 
makes the life-time slightly longer than that of the plain 
BSSN, while A-adjustment increases the life-time of the 
simulation by a factor of 10. However, it is also true that 
the error grows in time in all the three cases. 

We also find that the error is induced by distortion 
of the wave, i.e. the both phase and amplitude errors 
distort the numerical solution. In Fig. EKb), we show a 
snapshot of "f xx numerical solution at T = 100, together 
with the exact solution at the same time coordinate. The 
amplitude difference between the numerical and exact so- 
lutions is apparently less when we use the A-eq. adjusted 
system than that of the plain system. In Sec. I VII later, 
we discuss what causes the error and why the simula- 
tion life-time becomes longer when we use the adjusted 
system. 
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FIG. 1: The one-dimensional gauge-wave test with the plain BSSN system. The L2 norm of H and M x , rescaled by p 2 /4, are 
plotted with a function of the crossing-time. The amplitude of the wave is A = 0.01. The loss of convergence at the early time, 
near the 20 crossing-time, can be seen, and it will produce the blow-ups of the calculation in the end. 
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FIG. 2: The one-dimensional gauge- wave test with the adjusted BSSN system in the ^4-equation (|2.20[) . The L2 norm of Ti 
and M x , rescaled by p 2 /4, are plotted with a function of the crossing-time. The wave parameter is the same as with Fig. [TJ 
and the adjustment parameter ka is set to ha = 0.005. We see the higher resolution runs show convergence longer, i.e., the 
300 crossing-time in Ti and the 200 crossing-time in M x with p = 4 and 8 runs. All runs can stably evolve up to the 1000 
crossing-time. 



B. Linear wave test 



The second test is the linear wave propagation test, 
§111 B, to check the accuracy of wave propagations in the 
adjusted systems. We find that the linear wave testbed 



does not produce a significant constraint violation even 
for the plain BSSN system. The simulation does not 
crash at the 1000 crossing-time irrespective of the res- 
olutions. Figure [5] illustrates the profiles of ^ zz — 1 at 
the 500 crossing-time. The figure indicates the simula- 
tion does not produce the amplitude error but does pro- 
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FIG. 3: The one-dimensional gauge- wave test with the adjusted BSSN system in the T-equation p. 221) . The L2 norm of TL 
and M Xl rescaled by p 2 /4, are plotted with a function of the crossing-time. The wave parameter is the same as Fig.Q] and the 
adjustment parameter is Kf = —0.1. Note the near perfect overlap for the 55 crossing-time in TL and the 40 crossing-time in 

Mm. 
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FIG. 4: Evaluation of the accuracy of the one-dimensional gauge-wave testbed. Lines show the plain BSSN, the adjusted BSSN 
with „4-equation, and with F-equation. (a) The L2 norm of the error in *y xx , using (|4,2[) . (b) A snapshot of the exact and 
numerical solution at T = 100. 



duce the phase error. However, we also observe that the 
higher resolution run reduces the phase error. We tried 
the same evolutions with adjusted BSSN systems. How- 
ever, all the results are indistinguishable from the those 
of the plain BSSN system. This is because the adjusted 
terms of the equations are small due to the small viola- 
tions of constraints. Figure [6] shows a snapshot of the 
error defined by 7 Z2 — 7l° xact ^ at the 500 crossing-time 
both for the plain BSSN system and the adjusted BSSN 
system where the A-equation where ka = 10~ 3 . Since 



two lines are matching quite well, we can say that the 
adjusted BSSN system produces the same result as the 
plain BSSN system, including the phase error. Results 
from the other adjusted BSSN systems are almost the 
same qualitatively, including their convergence features. 
We also remark that we do not see a case in which ad- 
justment worsens accuracy and stability. 
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FIG. 5: Snapshots of the one-dimensional linear wave at different resolutions with the plain BSSN system at the simulation 
time 500 crossing-time. We see there exists phase error, but they are convergent away at higher resolution runs. 
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FIG. 6: Snapshot of errors with the exact solution for the Linear Wave testbed with the plain BSSN system and the adjusted 
BSSN system with the A equation at T — 500. The highest resolution p = 8 is used in both runs. The difference between the 
plain and the adjusted BSSN system with the A equation is indistinguishable. Note that the maximum amplitude is set to be 
10~~ 8 in this problem. 



C. Gowdy-wave test 

The third test is the polarized Gowdy-wave test, §111 
C, to check the adjustments in the strong field regime. 



The plain BSSN 



reported they can produce the stable and accurate evo- 
lution for the 1000 crossing-time by implementing the 
higher order differencing scheme to their LazEv code. 
However, it should be emphasized that they suggested 
their code produces the stable simulation only when they 
used the Kreiss-Oliger dissipation term [3(|. ) 



In Fig. [7J We first show the case of the plain BSSN 
evolution. We find that the second-order convergence 
continues up-to the 100 crossing-time and the higher res- 
olutions runs tend to crash at early times. This behavior 
(and crashing time) almost coincides with the results of 
the Cactus BSSN code, reported by Alcubierre et al. 29 1 
(see their Fig. 7). (We remark that Zlochower et al. 32 1 



2. Adjusted BSSN with A- equation 

Adjustment of the A-equation using the momentum 
constraint (|2.20[> . extends the life-time of the simulation 
10 times longer for the highest resolution run. Figure [5] 
depicts the rescaled L2 norm of Tt and M. z versus time. 
We set Kj\ = —0.001. (Note that the signature of k is 
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FIG. 7: Collapsing polarized Gowdy-wave test with the plain BSSN system. The L2 norm of H and Mz, rescaled by p /4, are 
plotted with a function of the crossing-time. (Simulation proceeds backwards from t = 0.) We see almost perfect overlap for 
the initial 100 crossing-time, and the higher resolution runs crash earlier. This result is quite similar to those achieved with 
the Cactus BSSN code, reported by [29| . 



reversed from the expected one, since the evolution is 
backward in time.) 

We find that an almost perfect overlap up to the 1000 
crossing-time under both the Hamiltonian constraint and 
the Momentum constraint. (These overlaps indicate that 
the error in H and M. z in the p = 8 resolution runs are 
16 times smaller than these errors in the p = 2 resolution 
run. ) However, we also find oscillations in the Momen- 
tum constraint, especially in the end of the simulation. 



3. Adjusted BSSN with "/-equation 

The case of the adjustment of the 7-equation using the 
(/-constraint (|2.2ip . is shown in Fig. [9] The adjustment 
parameter Kj is set at 0.000025. (Again, the signature of 
k is reversed from the expected one.) 

Figure [5] shows that an almost perfect overlap is ob- 
tained for the 200 crossing-time in both Tt and Ai z . 
The higher resolution runs tend to crash at earlier times, 
which is same as with the plain BSSN system. However, 
the convergence time becomes longer than that of the 
plain BSSN system. We will discuss the quantitative im- 
provement for the 7-adjustment in the next subsection. 



4- Adjustment effect 

In order to check the accuracy of the simulations, we 
prepare Fig. [TU] to show the error of the j zz component 



of the metric. 

Unlike the gauge- wave or the linear wave test, in this 
Gowdy-wave test the amplitude of the metric functions 
damps with time. Therefore we use the criterion that 
the error normalized by j zz be under 1% for an accurate 
evolution. This criterion is the same as the one used in 
Zlochower et al. [HJ ■ 

Figurc[Tn]shows the normalized error in "f zz versus time 
for the plain BSSN, adjusted BSSN with ^-equation, and 
adjusted BSSN with 7-equation systems. We find that 
these three systems produce accurate results up to t = 
200, t = 1000, and t = 400, respectively. This proves that 
the adjustments work effectively, i.e, they make possible a 
stable and accurate simulation, especially the .4-adjusted 
BSSN system. 



VI. SUMMARY AND DISCUSSION 

In this article, we presented our numerical compar- 
isons of the BSSN formulation and its adjusted versions 
using constraints. We performed three testbeds: gauge- 
wave, linear wave, and collapsing polarized Gowdy-wave 
tests with their evolutions by three kinds of adjust- 
ments, which were previously proposed by Yoneda and 
Shinkai [2(1 ] based on their constraint propagation anal- 
ysis. 

The idea of the adjusted systems is to construct a sys- 
tem robust against constraint violations by modifying the 
evolution equations using the constraint equations. 

We can summarize our tests as follows: 
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FIG. 8: Collapsing polarized Gowdy-wave test with the adjusted BSSN system in the j4-equation (|2.20|) . with ka = —0.001. 
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the initial 1000 crossing-time in both constraint equations, Ti and M z , even for the highest resolution run. 
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FIG. 9: Collapsing polarized Gowdy-wave test with the adjusted BSSN system in the 7-equation (|2.21|) . with = 0.000025. 
The figure style is the same as Figure [7J Note the almost perfect overlap for 200 crossing-time in the both the Hamiltonian 
and Momentum constraint and the p = 2 run can evolve stably for 1000 crossing-time. 



• When the plain (original) BSSN evolutions already 
show satisfactory good evolutions (e.g., the linear 
wave test), the constraint violations (i.e., adjusted 
terms) are also small or ignorable. 

Therefore the adjusted BSSN equations become quite 
similar to the plain BSSN equations, and their results 
coincide with the plain BSSN results. 



• Among the adjustments we tried, we observed that 
the adjusted BSSN system with the A-eq. (1^0]) 
is the most robust for all the testbeds examined in 
this study. It gives us an accurate and stable evo- 
lution compared to the plain BSSN system. Quan- 
titatively, the life-time of the simulation becomes 
10 times longer for the gauge- wave testbed and 5 
times longer for the Gowdy-wave testbed than the 



12 



o 




1e-08 
-1 



000 -800 -600 -400 -200 
t(harmonic) 



200 



FIG. 10: Comparisons of systems in the collapsing polarized Gowdy-wave test. The L2 norm of the error in 7 2Z , rescaled by the 
L2 norm of y zz , for the plain BSSN, adjusted BSSN with ^4-equation, and with 7-equation are shown. The highest resolution 
run, p = 8, is depicted for the plots. We can conclude that the adjustments make longer accurate runs available. Note that the 
evolution is backwards in time. 



life-time of the plain BSSN system. However, it 
should be noted that for the gauge- wave testbed, 
the convergence feature is lost at a comparatively 
early time, the 200 crossing-time in the Hamilto- 
nian constraint and the 50 crossing-time in the mo- 
mentum constraint. 

Recently, it has been claimed that the set up of the gauge 
wave pro blem in Apples-with-Apples has a problematic 
point [37l| . w hich arises from the harmonic gauge condi- 
tion. In [41(, it is argued that this gauge has a residual 
freedom in the form H — > e xt H, where A is an arbitrary 
and H is a function in Eq. (|3.1|) . Of course, our set up 
corresponds to the A = case, but numerical error easily 
excites modes that result in either exponentially increas- 
ing or decaying metric amplitude. Actually, we find the 
amplitude of the error decays with time in this testbed. 
So, we conclude that due to the adjustment, the growing 
rate of the gauge mode is suppressed and the life-time of 
the simulation is extended as a result. 

• The other type of adjustments (|2~2"T1 and |2~2"2"|) 

show their apparent effects while depending on a 
problem. The T-adjustment for the gauge-wave 
testbed makes the life-time longer slightly. The 7- 
adjustment for the Gowdy-wave testbed makes pos- 
sible a simulation twice as long as the plain BSSN 
system. 

We can understand the effect of the adjustments in terms 
of adding dissipative terms. By virtue of the definition 
of the constraints, we can recognize that the adjusted 
equation corresponds to the diffusion equation (see, for 
example, Eq. (|2.20p ) and the signature of k determines 
whether the diffusion is positive or negative. In the ad- 
justed A-eq. system, (|2.20l) . the adjustment term corre- 
sponds to the positive diffusive term, due to the defini- 



tion of Mi and the positiveness of ka (see Eq. (|2.15l) and 
(|2.20p ). This fact might explain why the adjusted A-eq. 
system works effectively for all the testbeds. 

In contrast, why are not all the adjustments effective 
in all testbeds? As we mentioned in Sec. IIB, the eigen- 
value analysis was made on the linearly perturbed viola- 
tion of constraints on the Minkowski space-time. Since 
the constraint violation grows non-linearly as seen in the 
Appendix of (2d| . the candidates may not be the best in 
their later evolution phase. 

We remark upon two more interesting aspects arising 
from our study. The first is the mechanism of the con- 
straint violations. As was shown in the appendix of [26| . 
each constraint propagation (behavior of their growth or 
decrease) depends on the other constraint terms together 
with itself. That is, we can guess A and S constraints 
(|2.17l and I2.18[) in this article, propagate independently 
of the other constraints, while the violation of the Q- 
constraint, (|2.16[) is triggered by the violation of the mo- 
mentum constraint, and both the Hamiltonian and the 
momentum constraints are affected by all the other con- 
straints. Such an order of the constraint violation can be 
guessed in Fig. [11] (earlier time), where we plot the rate 
of constraint violation normalized with its initial value, 
1^112(0/11^112(0), as a function of time, for the gauge- 
wave testbeds with the plain BSSN evolution. (Note that 
the constraints at the initial time, <5C(0), are not zero due 
to the numerical truncation error. ) The parameters are 
the same as those shown in Sec. MI A[ and the lowest 
resolution run is used. From this investigation, we might 
conclude that to monitor the momentum constraint vio- 
lation is the key to checking the stability of the evolution. 

The second remark is on the Lagrange multipliers, k, 
used in the adjusted systems. As discussed in Sec. Ill Bt 
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FIG. 11: The violation of all constraints normalized with their initial values, ||<5C||2(i)/||<5C||2(0), are plotted with a function 
of time. The evolutions of the gauge-wave testbeds with the plain BSSN system are shown. The parameters of the test are 
the same as those shown in Sec. MI Al and the lowest resolution run, p = 2, is applied. By observing which constraint triggers 
the other constraint's violation from the constraint propagation equations, we may guess the mechanism by which the entire 
system is violating accuracy and stability. See the text for details. 



the signatures of the ks are determined a priori, and we 
confirmed that all the predicted signatures of ks in [26[ 
are right to produce positive effects for controlling con- 
straint violations. However, we have to search for a suit- 
able magnitude of ks for each problem. Therefore we are 
now trying to develop a more sophisticated version, such 
as an auto-controlling n system, which will be reported 
upon in the future elsewhere. 

Although the testbeds used in this work are simple, it 
might be rather surprising to observe the expected effects 
of adjustments with such a slight change in the evolution 
equations. We therefore think that our demonstrations 
imply a potential to construct a robust system against 
constraint violations even in highly dynamical situations, 
such as black hole formation via gravitational collapse, 
or binary merger problems. We are now preparing our 
strong-field tests of the adjusted BSSN systems using 
large amplitude gravitational waves, black hole space- 



time, or non-vacuum space-time, which will be reported 
on in the near future. 
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We present our numerical comparisons between the BSSN formulation widely used in numerical 
relativity today and its adjusted versions using constraints. We performed three testbeds: gauge- 
wave, linear wave, and Gowdy-wave tests, proposed by the Mexico workshop on the formulation 
problem of the Einstein equations. We tried three kinds of adjustments, which were previously 
proposed from the analysis of the constraint propagation equations, and investigated how they 
improve the accuracy and stability of evolutions. We observed that the signature of the proposed 
Lagrange multipliers are always right and the adjustments improve the convergence and stability of 
the simulations. When the original BSSN system already shows satisfactory good evolutions (e.g., 
linear wave test), the adjusted versions also coincide with those evolutions; while in some cases (e.g., 
gauge-wave or Gowdy-wave tests) the simulations using the adjusted systems last 10 times as long 
as those using the original BSSN equations. Our demonstrations imply a potential to construct a 
robust evolution system against constraint violations even in highly dynamical situations. 

PACS numbers: 04.20.-q, 04.25. Dm, 04.25.-g 



I. INTRODUCTION 

Numerical integration of the Einstein equations is the 
only way to investigate highly dynamical and nonlinear 
gravitational space-time. The detection of gravitational 
wave requires templates of waveform, among them merg- 
ers of compact objects are the most plausible astrophysi- 
cal sources. Numerical relativity has been developed with 
this purpose over decades. 

For neutron star (NS) binaries, a number of scien- 
tific numerical simulations have been done so far, and 
we are now at the level of discussing the actual physics 
of the phenomena, including the effects of the equa- 
tions of state, hydrodynamics, andgeneral relativity by 
evolving various initial data [H, UHi H HI- Mergers 
of black holes (BHs) are also available after the break- 
through by Pretorius [6j in 2004. Pretorius's implemen- 
tation had many novel features in his code; among them 
he discretizes the four-dimensional Einstein equations di- 
rectly, which is not a conventional approach so far. How- 
ever, after the announcements of successful binary BH 
mergers by Campanelli et al. [3] and Baker et al. @ 
based on the standard 3+1 decomposition of the Einstein 
equations, many groups be gan producing interesting re- 
sults [I EE El El EE EllS ES E3, UM- The merger 
of NS-BH binary simulations has also been reported re- 
cently, e.g. 

Almost all the groups which apply the above conven- 
tional approach use the so-called BSSN variables together 
with "1 + log" -type slicing conditions for the lapse func- 
tion and 'T-driver" type slicing conditions for the shift 
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function. BSSN stands for Baumgarte-Shapiro [2fJ and 
Shibata-Nakamura [21| . the modified Arnowitt-Deser- 
Misner formulation initially proposed by Nakamura [22] • 
(The details are described in i jll Al ) There have al- 
ready been several efforts to explain why the combi- 
nation of this recipe works from the point of view of 
the well- pos edness of the partial differential equations 
(e.g. [23L l24j). However, the question remains whether 
there exists an alternative evolution system that enables 
more long-term stable and accurate simulations. The 
search for a better set of equations for numerical inte- 
grations is called the formulation problem for numerical 
relativity, of which earlier stages are reviewed by one of 
the authors [251 ]. 

In this article, we report our numerical tests of mod- 
ified versions of the BSSN system, the adjusted BSSN 
systems, proposed by Yoneda and Shinkai [26] . The idea 
of their modifications is to add constraints to the evolu- 
tion equations like Lagrange multipliers and to construct 
a robust evolution system which evolves to the constraint 
surface as the attractor. Their proposals are based on the 
eigenvalue analysis of the constraint propagation equa- 
tions (the evolution equations of the constraints) on the 
perturbed metric. For the ADM formulation, they ex- 
plain why the standard ADM does not work for long-term 
simulations by showing the existence of the constraint vi- 
olating mode in perturbed Schwarzschild space-time [27| . 
For the BSSN formulation, they analyzed the eigenval- 
ues of the constraint propagation equations only on flat 
space-time [261 ] . but one of their proposed adjustments 
was immediately tested by Yo et al. [28j for the numeri- 
cal evolution of Kerr-Schild space-time and confirmed to 
work as expected. (The details arc described in SQlBj) 

Our numerical examples are taken from the proposed 
problems for testing the formulations of the Mexico 
Numerical Relativity Workshop 2001 participants [29| . 
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which are sometimes called the Apples- with- Apples test. 
To concentrate the comparisons on the formulation prob- 
lem, the templated problems are settled so as not to 
require technical complications; e.g., periodic boundary 
conditions are used and the slicing conditions do not re- 
quire solving elliptical equations. Several groups already 
reported their code tests using these Apples tests (e.g. 
[3Q , I3II , H2| ) , and we are also able to compare our results 
with theirs. 

This article is organized as follows. We describe the 
BSSN equations and the adjusted BSSN equations in 
Sec. Ill Al and III Bl We give our three numerical test 
problems in Sec. IIIII Comments on our coding stuff are 
in Sec. [IVj Sec. |V] is devoted to showing numerical re- 
sults for each testbeds, and we summarize the results in 
Sec. ED 



widely used among numerical relativists. 

The idea of the BSSN formulation is to introduce aux- 
iliary variables to those of the Arnowitt-Deser-Misner 
(ADM) formulation for obtaining longer stable numerical 
simulations. The basic variables of the BSSN formulation 
are (0, 7^ , K, Aij , T*), which are defined by 



= Y2 log ( det7i ^' 



Hi 
K 



A 



— c -40 



Ki 



1 jki 



(2.1) 

(2.2) 
(2.3) 

(2.4) 
(2.5) 



II. BASIC EQUATIONS 

A. BSSN equations 

We start by presenting the standard BSSN formula- 
tion, where we follow the notations of [20(, which are 



where (7^ , Kij ) are the intrinsic and extrinsic ADM 3- 
metric. The conformal factor <p is introduced so as to 
set 7 = det[7y] as unity, Aij is supposed to be trace- 
less, and r l is treated independently in evolution equa- 
tions. Therefore these three requirements turn into the 
new constraints [below (|2T6|l - (|238|) ]. 

The set of the BSSN evolution equations are 



M = -\aK + /3 l d l( f>+ld t f3\ 
6 6 

2 

drfij = -laAij + % k d 3 /3 k + j jk di/3 k - -%d k p k + /3 k d k %, 
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(2.6) 
(2.7) 
(2.8) 

(2.9) 
(2.10) 



where Di is the covariant derivative with respect to the 3- metric 7y and TF means trace- free operation, i.e., Hj F 
Hij — ^jijH k k . The Ricci tensor is computed with the conformal connection T l as 



R% = -2A£»j0 - 2^D k D k 4> + 4D l( j>D^ - A%D k 4>D k ^ 
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(2.11) 
(2.12) 

(2.13) 



where Di is a covariant derivative associated with 7^ . the Hamiltonian and momentum constraint equations, 
Similarly to the ADM formulation, this system has are expressed in terms of the BSSN basic variables and 
constraint equations. The two "kinematic" constraints, are written as 
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Additionally, the BSSN formulation requires three "alge- 
braic" constraint relations; 
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r 

7-lwO, 



0, 



(2.16) 
(2.17) 
(2.18) 



where f|2 . and (|2 . 1 7[) are from the definitions of 
and (|2.4|) . respectively. Equation (|2.18|) is from the re- 
quirement on 7. 

These five constraints are, theoretically, supposed to be 
zero at all times; therefore they can be used to modify the 
dynamical equations. For example, Alcubierre et al. [34| 
announced that the replacement of the terms in (|2.10j) 
using the momentum constraint drastically changes the 
stability feature. Actually, such replacements of terms 
using constraints are applied (with/without intentions) 
in many terms in (|2.6D - (|2.10[) . which are expressed as 
Eqs. (2.27)-(2.31) in j2j|. 

Alcubierre et al. [35| also pointed out that the re- 
definition of Aij by 



A, 



A,. 



trA 



(2.19) 



during the time evolutions improves the numerical sta- 
bility. This technique again can be understood as the 
trace-out of the .A-constraint (|2.17|) from the evolution 
equations. In our numerical code, we do not apply this 
technique because we recognize the trace-free property 
as the new constraint A in the BSSN system, and our 
purpose is to construct a system preventing the violation 
of constraints. 

Recently, several groups applied artificial dissipation 
(e.g. [36] ]) to obtain stable evolutions (see, e.g. [HI, [H, 
137]). We, however, do not introduce such dissipations 
in our code, since we try to clarify the difference of sta- 
bility from the viewpoint of formulations of the Einstein 
equations. 



B. Adjusted BSSN systems 

To understand the stability property of the BSSN sys- 
tem, Yoneda and Shinkai [2(| studied the structure of the 
evolution equations, (|2.6p - (|2.10p , in detail, especially how 
the modifications using the constraints, (|2.14p - (|2.18p . af- 
fect to the stability. They investigated the signature of 
the eigenvalues of the constraint propagation equations 
(dynamical equations of constraints), and explained that 



the standard BSSN dynamical equations are balanced 
from the viewpoints of constrained propagations, includ- 
ing a clarification of the effect of the replacement using 
the momentum constraint equation. 

Moreover, they predicted that several combinations 
of modifications have a constraint-damping nature, and 
named them adjusted BSSN systems. (Their predictions 
are based on the signature of eigenvalues of the constraint 
propagations, and the negative signature implies a dy- 
namical system which evolves toward the constraint sur- 
face as the attractor.) 

Among them, in this work, we test the following three 
adjustments: 

1. An adjustment of the A-equation with the momen- 
tum constraint: 



d t A 



K A aD(iMj), 



(2.20) 



where ka is predicted (from the eigenvalue analy- 
sis) to be positive in order to damp the constraint 
violations. 

2. An adjustment of the 7-equation with Q constraint: 
daij = + ^a% (i D j) g k , (2.21) 



with k~ < 0. 



3. An adjustment of the T-equation with Q constraint: 
d t f l = dft i + K t aQ i . (2.22) 



with Kp < 0. 



These three adjustments are chosen as samples of "best 
candidates", Eq. (4.9)-(4.11) in 26]. The term "best" 
comes from their conjecture on the eigenvalue analysis 
of the constraint propagation matrix; that is, (a) all the 
resultant eigenvalues from above adjustments can be less 
than or at most equal to zero, which indicates the de- 
cay of constraint errors, and (b) the resultant constraint 
propagation matrix is diagonalizable, which guarantees 
the predictions of above eigenvalue analysis (see Table II 
in [26j). However, since above eigenvalues include zero 
elements and also above analysis assumes a linearly per- 
turbed metric about the flat space-time, the effects of 
the adjustments (|2.20p - (|2.22p need to be demonstrated 
via numerical experiments. 
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III. NUMERICAL TESTBED MODELS 



where 



Following the proposals of the Mexico Numerical Rela- 
tivity Workshop [29j, we perform three kinds of tests. In 
this section, we explicitly give some details of the models. 



A. Gauge-wave testbed 

The first test is the trivial Minkowski space-time, but 
sliced with the time-dependent 3-metric, which is called 
the gauge- wave test. The 4- metric is obtained by coordi- 
nate transformation from the Minkowski metric as 



ds 2 



-Hdt 2 + Hdx 2 + dy 2 + dz 2 , 



where 



H = H(x -t) = l- Asin 



2tt(x - t) 



(3.1) 



(3.2) 



which describes a sinusoidal gauge wave of amplitude A 
propagating along the x-axis. The non-trivial extrinsic 
curvature is 



K, 



irA 



1 + A sin 



2ir(x-t) 



(3.3) 



Following [29l | . we chose numerical domain and parame- 
ters as follows: 

• Gauge- wave parameters: d — 1 and A = 10~ 2 

• Simulation domain: 0.5, 0.5], y = z = 

• Grid: x l = —0.5 + (n — ^)dx with n = 1, • • • 50p, 
where dx = l/(50p) with p = 2,4, 8 

• Time step: dt = 0.25dx 

• Boundary conditions: Periodic boundary condition 
in x direction and planar symmetry in y and z di- 
rections 



• Gauge conditions: 

d t a = -a 2 K, ft = 0. 



(3.4) 



The ID simulation is carried out for a T — 1000 crossing- 
time or until the code crashes, where one crossing-time 
is defined by the length of the simulation domain. 



B. Linear wave testbed 

The second test is to check the ability of handling a 
travelling gravitational wave. The initial 3-metric and 
extrinsic curvature Kij are given by a diagonal pertur- 
bation with component 



b = A sin 



2tt(x - t) 



(3.6) 



for a linearized plane wave traveling in the ir-direction. 
Here d is the linear size of the propagation domain and A 
is the amplitude of the wave. The non-trivial components 
of extrinsic curvature arc then 



1 



1 



K yy = ~7;dtb, K zz = -d t b. 



(3.7) 



2 2 
Following [2^|, we chose the following parameters: 

• Linear wave parameters: d = 1 and A = 10~ 8 

• Simulation domain: xg[~ 0.5, 0.5], y — 0, z = 

• Grid: x l = —0.5 + (n — ^)dx with n = 1, • • • 50p, 
where dx = l/(50p) with p = 2, 4, 8 

• Time step: dt = 0.25dx 

• Boundary conditions: Periodic boundary condition 
in x direction and planar symmetry in y and z di- 
rections 

• Gauge conditions: a = 1 and (3 l — 

The ID simulation is carried out for a T = 1000 crossing- 
time or until the code crashes. 



C. Collapsing polarized Gowdy-wave testbed 

The third test is to check the formulation in a strong 
field context using the polarized Gowdy metric, which is 
written as 

ds 2 = t- l ' 2 e x ' 2 {-dt 2 + dz 2 ) + t(e p dx 2 + e- p dy 2 )(3.8) 

Here time coordinate t is chosen such that time increases 
as the universe expands. Simple forms of the solutions, 
P and A, are given by 

P = J (27rf)cos(27rz), (3.9) 
A = -27rfJ (27rf)Ji(27rf)cos 2 (27rz) 

+ 27T 2 t 2 [J 2 (27rf) + J 2 (27rf)] 

1, 



^ (27r) 2 [J 2 (27r) + J 2 (2^)] 
-2^(2^)^(2^)], 



(3.10) 



where J n is the Bessel function. The non-trivial extrinsic 
curvatures are then 



ds 2 



-dt 2 + dx 2 + (1 + b)dy 2 + (1 - b)dz 2 , (3.5) 



K 



K 7 . 



_l t V4 e -A/4 e P (1+tP)t)) (3 . n) 

_ItV4 e - A /4 e --P(i_tP il ) i (3.12) 

It-V4 e A/4 (t -l_ At) (313) 
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According to [29j , the new time coordinate t, which satis- 
fies harmonic condition, is obtained by coordinate trans- 
formation as 

t(r) = ke CT , (3.14) 

where c and k are arbitrary constants. Using this free- 
dom, we can set the lapse function in the new time co- 
ordinate to be unity at the initial time. Concretely, we 
set 

to = r ~ 9.8753205829098, 

c - 0.0021195119214617, (3.15) 
k ~ 9.6707698127638, 

where to is the initial time. Following [29l |. we perform 
our evolution in the collapsing (i.e. backward in time) 
direction. Parameters are chosen as follows: 

• Simulation domain: z £ [—0.5, 0.5], x = y = 

• Grid: z — —0.5 + (n — ^)dz with n = 1, • • • 50p, 
where dz = l/(50p) with p = 2,4,8 

• Time step: dt = 0.25dz 

• Boundary conditions: Periodic boundary condition 
in z-direction and plane symmetry in x- and y- 
directions 

• Gauge conditions: the harmonic slicing (|3.4j) and 

/& = 

The ID simulation is carried out for a T — 1000 crossing- 
time or until the code crashes. 



IV. THE CODE 
A. Code description 

We have developed a new numerical code based 
on the adjusted BSSN systems. The variables are 
(</>, Jij,K, Aij, T l ), and the evolution equations are (|2.6p ~ 
(f2~TUj) with/without adjustment (f2T20|) . (f2~2Tj) , and/or 
(|2.22[) . The time-integration is under the free- evolution 
scheme, and we monitor five constraints, (|2.14|) - (|2.18|) , to 
check the accuracy and stability of the evolutions. 

Our time-integration scheme is the three-step iterative 
Crank-Nicholson method with centered finite difference 
in space [39l | . This scheme should have second-order con- 
vergence both in space and time, and we checked its con- 
vergence in all the testbeds. 

As we have already mentioned in the end of §11 A, we 
do not apply the trace-out technique of , (|2.19p in our 
code. 

We also remark on our treatment of the conformal con- 
nection variable f 1 . As was pointed out in [H| , it is better 
not to use T l in all the evolution equations. We surmise 
this is because the amplification of the error due to the 



discrepancy of the definition (|2.5p . i.e., the accumula- 
tions of the violations of ^'-constraint (|2.16[) . Therefore, 
we used the evolved T 1 only for the terms in (|2.10|) and 
(|2.13[) . and not for other terms, so as not to implicitly 
apply the ^'-constraint in time evolutions. 



B. Debugging procedures 

It is crucial that our code can produce accurate re- 
sults, because the adjustment methods are based on 
the assumption that the code represents the BSSN sys- 
tem (|2.6p - (|2.10p accurately. We verified our code by com- 
paring our numerical data with analytic solutions from 
the gauge-wave and Gowdy-wave testbeds in Sec. IIIII 
The actual procedures are as follows: 

1. Evolve only one component, e.g. A xx , numerically, 
and express all the other components with those of 
the analytic solution. In this situation, the origin 
of the error is from the finite differencing of the 
analytic solution in the spatial direction and from 
that of the numerically evolved component (A xx ) 
both in spatial and time directions. We checked 
the code by monitoring the difference between the 
numerically evolved component (A xx ) and its ana- 
lytic expression. This procedure was applied to all 
the components one by one. 

2. Evolve only several components, e.g., A xx and T x , 
numerically, and express the other components by 
the analytic solution. The error can be checked by 
a procedure similar to the one above. 

3. Evolve all the components numerically, and check 
the error with the analytic solution. 

We repeated these procedure three times by switch- 
ing the propagation directions (x, y, and z-directions) 
of gauge-wave and Gowdy-wave solutions. We also ap- 
plied these procedures in a 2D test [29[ , and checked the 
off-diagonal component. 



C. Error evaluation methods 

It should be emphasized that the adjustment effect has 
two meanings, improvement of stability and of accuracy. 
Even if a simulation is stable, it does not imply that the 
result is accurate. We judge the stability of the evolution 
by monitoring the L2 norm of each constraint, 



\\SCMt)= /-£(C(t; W )) 2 , (4.1) 

where N is the total number of grid points, while 
we judge the accuracy by the difference of the met- 
ric components gij (t; x, y, z) from the exact solution 
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7iJ xact ' ) (i; x, y, z), higher resolution cases. However, it is also true that all 

errors are still under the errors of the plain BSSN system. 
/ i / (exact) \ 2 Therefore, we conclude that this adjusted system shows 

\\^9ij\\2{t) = . — \ 9ij ~ 9ij J ■ (4-2) a weaker instability than the plain system. 



D. Magnitude of k 

Adjusted systems, (|2.20l) - (|2.22l) . require to specify the 
parameter n. From the analytical prediction in [261 ] we 
know the signature of k, but not for its magnitude. By 
definition of the adjustment terms in Eq. (|2.20[) - (|2.22[) . 
applying small k should produce the close results with 
those of the plain system. On the contrary, the large 
k system will violate the Courant-Friedrich-Lewy condi- 
tion (4^|. Hence, there exists a suitable region in the 
adjustment parameters. 

At this moment, we have to chose k experimentally, 
by observing the life-time of simulations. The value of /?, 
used in our demonstrations, is one of the choices of which 
the adjustment works effectively in all the resolutions. 

V. NUMERICAL RESULTS 
A. Gauge-wave test 

1. The plain BSSN system 

As the first test, we show the plain BSSN evolution 
(that is, no adjustments) in Fig. [T] for the gauge-wave 
test. In Fig. [TJ the L2 norms of the Hamiltonian and 
momentum constraints (|4. 1[) are plotted as a function of 
the crossing-time. The second-order convergent nature 
is lost at an early time, the 20 crossing-time, and the 
simulation crashes at about the 100 crossing-time. The 
poor performance of the plain BSSN system for the gauge 
wave test has been reported in [3l[ (see their Fig. 8). This 
drawback, on the other hand, can be overcome if one uses 
the fourth-order finite differencing scheme, an example of 
which can be seen in [33] (see their Fig. 2). 

2. Adjusted BSSN with A-equation 

We found that the simulation lasts 10 times longer 
with the adjustment in the A-equation using the mo- 
mentum constraint (|2.20[) . Figure [3] shows the L2 norms 
of the Hamiltonian and Momentum constraints in the 
same style as in Fig. [T] The adjustment parameter is set 
at ka = 0.005 for this plot. We obtain almost prefect 
overlap of the rescaled Hamiltonian constraint for 200 
crossing-times and almost perfect overlap in the momen- 
tum constraint for 50 crossing-times; there apparently 
improve the results of the plain BSSN system (see Fig.[T]). 
We show the plots until the 1000 crossing-time, there we 
observe the growth of the error both in later time and in 



3. Adjusted BSSN with T -equation 

The case of the adjustment of the T-equation using the 
Q constraint (|2.22[) is shown in Fig. [3] 

The adjustment parameter is set at Kf — —0.1. We 
find that the second-order convergence breaks down near 
the 40 crossing-time under the Momentum constraint, 
which is almost the same as with the plain BSSN system. 
However, the convergence of the Hamiltonian constraint 
is improved, i.e., it continues to the near 55 crossing-time. 
The life-time of the simulation is almost the same as that 
of the plain BSSN system. 



4- Adjusted BSSN with ^/-equation 

We also tested the cases of the adjustment of the 7- 
equation using the Q constraint, (|2.21[) . We again ob- 
served the effects of the adjustment on its stability and 
accuracy but found a rather small effect compared to the 
cases of the adjustments of (|2.20[) or (|2.22p . up to our 
trials of the parameter range of k 7 . Therefore we omit 
showing the results. 



5. Evaluation of Accuracy 

For evaluating the accuracy, we prepare Fig. Ufa), in 
which we plot the L2 norm of the error in j xx , (|4.2p , with 
the function of time. Three lines correspond to the result 
of the plain BSSN system, A-eq. adjusted, and f-eq. 
adjusted BSSN system, respectively. The T-adjustment 
makes the life-time slightly longer than that of the plain 
BSSN, while A-adjustment increases the life-time of the 
simulation by a factor of 10. However, it is also true that 
the error grows in time in all the three cases. 

We also find that the error is induced by distortion 
of the wave, i.e. the both phase and amplitude errors 
distort the numerical solution. In Fig. EKb), we show a 
snapshot of "f xx numerical solution at T = 100, together 
with the exact solution at the same time coordinate. The 
amplitude difference between the numerical and exact so- 
lutions is apparently less when we use the A-eq. adjusted 
system than that of the plain system. In Sec. I VII later, 
we discuss what causes the error and why the simula- 
tion life-time becomes longer when we use the adjusted 
system. 
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FIG. 1: The one-dimensional gauge-wave test with the plain BSSN system. The L2 norm of H and M x , rescaled by p 2 /4, are 
plotted with a function of the crossing-time. The amplitude of the wave is A = 0.01. The loss of convergence at the early time, 
near the 20 crossing-time, can be seen, and it will produce the blow-ups of the calculation in the end. 
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FIG. 2: The one-dimensional gauge- wave test with the adjusted BSSN system in the ^4-equation (|2.20[) . The L2 norm of Ti 
and M x , rescaled by p 2 /4, are plotted with a function of the crossing-time. The wave parameter is the same as with Fig. [TJ 
and the adjustment parameter ka is set to ha = 0.005. We see the higher resolution runs show convergence longer, i.e., the 
300 crossing-time in Ti and the 200 crossing-time in M x with p = 4 and 8 runs. All runs can stably evolve up to the 1000 
crossing-time. 



B. Linear wave test 



The second test is the linear wave propagation test, 
§111 B, to check the accuracy of wave propagations in the 
adjusted systems. We find that the linear wave testbed 



does not produce a significant constraint violation even 
for the plain BSSN system. The simulation does not 
crash at the 1000 crossing-time irrespective of the res- 
olutions. Figure [5] illustrates the profiles of ^ zz — 1 at 
the 500 crossing-time. The figure indicates the simula- 
tion does not produce the amplitude error but does pro- 
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FIG. 3: The one-dimensional gauge- wave test with the adjusted BSSN system in the T-equation p. 221) . The L2 norm of TL 
and M Xl rescaled by p 2 /4, are plotted with a function of the crossing-time. The wave parameter is the same as Fig.Q] and the 
adjustment parameter is Kf = —0.1. Note the near perfect overlap for the 55 crossing-time in TL and the 40 crossing-time in 
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FIG. 4: Evaluation of the accuracy of the one-dimensional gauge-wave testbed. Lines show the plain BSSN, the adjusted BSSN 
with „4-equation, and with F-equation. (a) The L2 norm of the error in *y xx , using (|4,2[) . (b) A snapshot of the exact and 
numerical solution at T = 100. 



duce the phase error. However, we also observe that the 
higher resolution run reduces the phase error. We tried 
the same evolutions with adjusted BSSN systems. How- 
ever, all the results are indistinguishable from the those 
of the plain BSSN system. This is because the adjusted 
terms of the equations are small due to the small viola- 
tions of constraints. Figure [6] shows a snapshot of the 
error defined by 7 Z2 — 7l° xact ^ at the 500 crossing-time 
both for the plain BSSN system and the adjusted BSSN 
system where the A-equation where ka = 10~ 3 . Since 



two lines are matching quite well, we can say that the 
adjusted BSSN system produces the same result as the 
plain BSSN system, including the phase error. Results 
from the other adjusted BSSN systems are almost the 
same qualitatively, including their convergence features. 
We also remark that we do not see a case in which ad- 
justment worsens accuracy and stability. 
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FIG. 5: Snapshots of the one-dimensional linear wave at different resolutions with the plain BSSN system at the simulation 
time 500 crossing-time. We see there exists phase error, but they are convergent away at higher resolution runs. 
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FIG. 6: Snapshot of errors with the exact solution for the Linear Wave testbed with the plain BSSN system and the adjusted 
BSSN system with the A equation at T — 500. The highest resolution p = 8 is used in both runs. The difference between the 
plain and the adjusted BSSN system with the A equation is indistinguishable. Note that the maximum amplitude is set to be 
10~~ 8 in this problem. 



C. Gowdy-wave test 

The third test is the polarized Gowdy-wave test, §111 
C, to check the adjustments in the strong field regime. 



The plain BSSN 



reported they can produce the stable and accurate evo- 
lution for the 1000 crossing-time by implementing the 
higher order differencing scheme to their LazEv code. 
However, it should be emphasized that they suggested 
their code produces the stable simulation only when they 
used the Kreiss-Oliger dissipation term [3(|. ) 



In Fig. [7J We first show the case of the plain BSSN 
evolution. We find that the second-order convergence 
continues up-to the 100 crossing-time and the higher res- 
olutions runs tend to crash at early times. This behavior 
(and crashing time) almost coincides with the results of 
the Cactus BSSN code, reported by Alcubierre et al. 29 1 
(see their Fig. 7). (We remark that Zlochower et al. 32 1 



2. Adjusted BSSN with A- equation 

Adjustment of the A-equation using the momentum 
constraint (|2.20[> . extends the life-time of the simulation 
10 times longer for the highest resolution run. Figure [5] 
depicts the rescaled L2 norm of Tt and M. z versus time. 
We set Kj\ = —0.001. (Note that the signature of k is 
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FIG. 7: Collapsing polarized Gowdy-wave test with the plain BSSN system. The L2 norm of H and Mz, rescaled by p /4, are 
plotted with a function of the crossing-time. (Simulation proceeds backwards from t = 0.) We see almost perfect overlap for 
the initial 100 crossing-time, and the higher resolution runs crash earlier. This result is quite similar to those achieved with 
the Cactus BSSN code, reported by [29| . 



reversed from the expected one, since the evolution is 
backward in time.) 

We find that an almost perfect overlap up to the 1000 
crossing-time under both the Hamiltonian constraint and 
the Momentum constraint. (These overlaps indicate that 
the error in H and M. z in the p = 8 resolution runs are 
16 times smaller than these errors in the p = 2 resolution 
run. ) However, we also find oscillations in the Momen- 
tum constraint, especially in the end of the simulation. 



3. Adjusted BSSN with "/-equation 

The case of the adjustment of the 7-equation using the 
(/-constraint (|2.2ip . is shown in Fig. [9] The adjustment 
parameter Kj is set at 0.000025. (Again, the signature of 
k is reversed from the expected one.) 

Figure [5] shows that an almost perfect overlap is ob- 
tained for the 200 crossing-time in both Tt and Ai z . 
The higher resolution runs tend to crash at earlier times, 
which is same as with the plain BSSN system. However, 
the convergence time becomes longer than that of the 
plain BSSN system. We will discuss the quantitative im- 
provement for the 7-adjustment in the next subsection. 



4- Adjustment effect 

In order to check the accuracy of the simulations, we 
prepare Fig. [TU] to show the error of the j zz component 



of the metric. 

Unlike the gauge- wave or the linear wave test, in this 
Gowdy-wave test the amplitude of the metric functions 
damps with time. Therefore we use the criterion that 
the error normalized by j zz be under 1% for an accurate 
evolution. This criterion is the same as the one used in 
Zlochower et al. [HJ ■ 

Figurc[Tn]shows the normalized error in "f zz versus time 
for the plain BSSN, adjusted BSSN with ^-equation, and 
adjusted BSSN with 7-equation systems. We find that 
these three systems produce accurate results up to t = 
200, t = 1000, and t = 400, respectively. This proves that 
the adjustments work effectively, i.e, they make possible a 
stable and accurate simulation, especially the .4-adjusted 
BSSN system. 



VI. SUMMARY AND DISCUSSION 

In this article, we presented our numerical compar- 
isons of the BSSN formulation and its adjusted versions 
using constraints. We performed three testbeds: gauge- 
wave, linear wave, and collapsing polarized Gowdy-wave 
tests with their evolutions by three kinds of adjust- 
ments, which were previously proposed by Yoneda and 
Shinkai [2(1 ] based on their constraint propagation anal- 
ysis. 

The idea of the adjusted systems is to construct a sys- 
tem robust against constraint violations by modifying the 
evolution equations using the constraint equations. 

We can summarize our tests as follows: 
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FIG. 8: Collapsing polarized Gowdy-wave test with the adjusted BSSN system in the j4-equation (|2.20|) . with ka = —0.001. 
The style is the same as in Fig. [7] and note that both constraints are normalized by p 2 /4. We see almost perfect overlap for 
the initial 1000 crossing-time in both constraint equations, Ti and M z , even for the highest resolution run. 
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FIG. 9: Collapsing polarized Gowdy-wave test with the adjusted BSSN system in the 7-equation (|2.21|) . with = 0.000025. 
The figure style is the same as Figure [7J Note the almost perfect overlap for 200 crossing-time in the both the Hamiltonian 
and Momentum constraint and the p = 2 run can evolve stably for 1000 crossing-time. 



• When the plain (original) BSSN evolutions already 
show satisfactory good evolutions (e.g., the linear 
wave test), the constraint violations (i.e., adjusted 
terms) are also small or ignorable. 

Therefore the adjusted BSSN equations become quite 
similar to the plain BSSN equations, and their results 
coincide with the plain BSSN results. 



• Among the adjustments we tried, we observed that 
the adjusted BSSN system with the A-eq. (1^0]) 
is the most robust for all the testbeds examined in 
this study. It gives us an accurate and stable evo- 
lution compared to the plain BSSN system. Quan- 
titatively, the life-time of the simulation becomes 
10 times longer for the gauge- wave testbed and 5 
times longer for the Gowdy-wave testbed than the 
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FIG. 10: Comparisons of systems in the collapsing polarized Gowdy-wave test. The L2 norm of the error in 7 2Z , rescaled by the 
L2 norm of y zz , for the plain BSSN, adjusted BSSN with ^4-equation, and with 7-equation are shown. The highest resolution 
run, p = 8, is depicted for the plots. We can conclude that the adjustments make longer accurate runs available. Note that the 
evolution is backwards in time. 



life-time of the plain BSSN system. However, it 
should be noted that for the gauge- wave testbed, 
the convergence feature is lost at a comparatively 
early time, the 200 crossing-time in the Hamilto- 
nian constraint and the 50 crossing-time in the mo- 
mentum constraint. 

Recently, it has been claimed that the set up of the gauge 
wave pro blem in Apples-with-Apples has a problematic 
point [37l| . w hich arises from the harmonic gauge condi- 
tion. In [41(, it is argued that this gauge has a residual 
freedom in the form H — > e xt H, where A is an arbitrary 
and H is a function in Eq. (|3.1|) . Of course, our set up 
corresponds to the A = case, but numerical error easily 
excites modes that result in either exponentially increas- 
ing or decaying metric amplitude. Actually, we find the 
amplitude of the error decays with time in this testbed. 
So, we conclude that due to the adjustment, the growing 
rate of the gauge mode is suppressed and the life-time of 
the simulation is extended as a result. 

• The other type of adjustments (|2~2"T1 and |2~2"2"|) 

show their apparent effects while depending on a 
problem. The T-adjustment for the gauge-wave 
testbed makes the life-time longer slightly. The 7- 
adjustment for the Gowdy-wave testbed makes pos- 
sible a simulation twice as long as the plain BSSN 
system. 

We can understand the effect of the adjustments in terms 
of adding dissipative terms. By virtue of the definition 
of the constraints, we can recognize that the adjusted 
equation corresponds to the diffusion equation (see, for 
example, Eq. (|2.20p ) and the signature of k determines 
whether the diffusion is positive or negative. In the ad- 
justed A-eq. system, (|2.20l) . the adjustment term corre- 
sponds to the positive diffusive term, due to the defini- 



tion of Mi and the positiveness of ka (see Eq. (|2.15l) and 
(|2.20p ). This fact might explain why the adjusted A-eq. 
system works effectively for all the testbeds. 

In contrast, why are not all the adjustments effective 
in all testbeds? As we mentioned in Sec. IIB, the eigen- 
value analysis was made on the linearly perturbed viola- 
tion of constraints on the Minkowski space-time. Since 
the constraint violation grows non-linearly as seen in the 
Appendix of (2d| . the candidates may not be the best in 
their later evolution phase. 

We remark upon two more interesting aspects arising 
from our study. The first is the mechanism of the con- 
straint violations. As was shown in the appendix of [26| . 
each constraint propagation (behavior of their growth or 
decrease) depends on the other constraint terms together 
with itself. That is, we can guess A and S constraints 
(|2.17l and I2.18[) in this article, propagate independently 
of the other constraints, while the violation of the Q- 
constraint, (|2.16[) is triggered by the violation of the mo- 
mentum constraint, and both the Hamiltonian and the 
momentum constraints are affected by all the other con- 
straints. Such an order of the constraint violation can be 
guessed in Fig. [11] (earlier time), where we plot the rate 
of constraint violation normalized with its initial value, 
1^112(0/11^112(0), as a function of time, for the gauge- 
wave testbeds with the plain BSSN evolution. (Note that 
the constraints at the initial time, <5C(0), are not zero due 
to the numerical truncation error. ) The parameters are 
the same as those shown in Sec. MI A[ and the lowest 
resolution run is used. From this investigation, we might 
conclude that to monitor the momentum constraint vio- 
lation is the key to checking the stability of the evolution. 

The second remark is on the Lagrange multipliers, k, 
used in the adjusted systems. As discussed in Sec. Ill Bt 
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FIG. 11: The violation of all constraints normalized with their initial values, ||<5C||2(i)/||<5C||2(0), are plotted with a function 
of time. The evolutions of the gauge-wave testbeds with the plain BSSN system are shown. The parameters of the test are 
the same as those shown in Sec. MI Al and the lowest resolution run, p = 2, is applied. By observing which constraint triggers 
the other constraint's violation from the constraint propagation equations, we may guess the mechanism by which the entire 
system is violating accuracy and stability. See the text for details. 



the signatures of the ks are determined a priori, and we 
confirmed that all the predicted signatures of ks in [26[ 
are right to produce positive effects for controlling con- 
straint violations. However, we have to search for a suit- 
able magnitude of ks for each problem. Therefore we are 
now trying to develop a more sophisticated version, such 
as an auto-controlling n system, which will be reported 
upon in the future elsewhere. 

Although the testbeds used in this work are simple, it 
might be rather surprising to observe the expected effects 
of adjustments with such a slight change in the evolution 
equations. We therefore think that our demonstrations 
imply a potential to construct a robust system against 
constraint violations even in highly dynamical situations, 
such as black hole formation via gravitational collapse, 
or binary merger problems. We are now preparing our 
strong-field tests of the adjusted BSSN systems using 
large amplitude gravitational waves, black hole space- 



time, or non-vacuum space-time, which will be reported 
on in the near future. 
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